In this lab you will begin to get used to using R with spatial data and practice your map making in R for a variety of data types.
Start by reading through sections 7.1-7.3 in Arnold and Tilton.2015. Humanities Data in R.
Then read through ‘Spatial Data Manipulation’ on RSpatial. [http://rspatial.org/spatial/index.html]. There’s a lot here, so feel free to move through it quickly and think of this as reference.
Upload to smartsite an Rmarkdown and HTML file.
Due: Monday, April 24, 9a
library (OpenStreetMap)
library (rgdal)
## Loading required package: sp
## rgdal: version: 1.2-5, (SVN revision 648)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
## Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
## Linking to sp version: 1.2-4
datapath = "/users/miadawson/Documents/~UC Davis~/~GEO 200CN/labs/lab 4" #assign directory location to an object called datapath
UCplaces <- read.csv(file.path(datapath, "UCplaces.csv")) #read UC places CSV
map <- openmap (c(38.55, -121.77), c(38.53, -121.74)) #provides the map boundary using upper-left and lower-right corners
mapLatLong <- openproj(map) #converts from the default mercator projection to Latitude Longitude
plot (mapLatLong)
points (UCplaces$long, UCplaces$lat, pch=16) #add the points
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3) #add labels
Create two additional maps: b) change the background context map, [hint: ?openmap will give you a list of possible types]
map2 <- openmap (c(38.55, -121.77), c(38.53, -121.74), type="stamen-toner") #change map type
mapLatLong2 <- openproj(map2)
plot (mapLatLong2)
points (UCplaces$long, UCplaces$lat, pch=16)
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3)
map3 <- openmap (c(38.57, -121.79), c(38.51, -121.71)) #adjust corners to zoom out on map
mapLatLong3 <- openproj(map3)
plot (mapLatLong3)
points (UCplaces$long, UCplaces$lat, pch=5, cex=.5) #change size and shape of points
text (UCplaces$long, UCplaces$lat, labels = UCplaces$names, pos=3, cex=.8) #change size of labels
library(sp)
library(maptools)
## Checking rgeos availability: TRUE
library(rgeos)
## rgeos version: 0.3-23, (SVN revision 546)
## GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921
## Linking to sp version: 1.2-4
## Polygon checking: TRUE
state <- readShapeSpatial(fn="State_2010Census_DP1") #read US Census data
## Warning: use rgdal::readOGR or sf::st_read
## Warning: use rgdal::readOGR or sf::st_read
index <- (as.data.frame(state)$STUSPS10 %in% c("AK", "HI")) #identify AK and HI data
state <- state[!index,] #remove AK and HI data
plot(state)
#convert projection to UTM
projectionObj <- CRS(projargs="+proj=longlat")
proj4string(state) <- projectionObj
stateTrans <- spTransform(x=state, CRSobj=CRS("+proj=utm +zone=14"))
#define centroid
centroid <- gCentroid(spgeom=stateTrans, byid=TRUE)
head(centroid)
## class : SpatialPoints
## features : 6
## extent : -1333889, 4104598, 2361966, 5482413 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=utm +zone=14 +ellps=WGS84
#plot map with labels placed at centroid
plot(stateTrans)
text(x=centroid$x, y=centroid$y, label=as.data.frame(stateTrans)$NAME10, cex=0.7)
and c) a cross hatch map showing the proportion of seasonal of recreational housing in each state.
#calculate proportion of housing units (DP0180001) declared as being “seasonal, recreational, or occasional” (DP0180008)
stateTransData <- as.data.frame(state)
perHouseRec <- stateTransData$DP0180008 / stateTransData$DP0180001
#convert to quantized IDs
bins <- quantile(perHouseRec, seq(0,1,length.out=8))
binId <- findInterval(perHouseRec, bins)
densityVals <- seq_len(length(bins)) * 5
plot(stateTrans, density=densityVals[binId])
You can also create a more familiar choropleth map with solid colors by defining your own color ramp/palette, or using one of the ramps defined within R.
colSet01 = c("grey90", "grey80", "grey70", "grey60", "grey50", "grey40", "grey30", "grey20")
plot (stateTrans, col=colSet01[binId])
colSet02 = rev(heat.colors(length(bins)))
plot (stateTrans, col=colSet02[binId])
#calculate proportion of population (DP0080001) declared as being “Black or African American alone or in combination with one or more other races ” (DP0090002)
pctBlack <- stateTransData$DP0090002 / stateTransData$DP0080001
#convert to quantized IDs
bins2 <- quantile(pctBlack, seq(0,1,length.out=8))
binId2 <- findInterval(pctBlack, bins2)
densityVals2 <- seq_len(length(bins2)) * 5
plot (stateTrans, col=colSet01[binId2])
plot (stateTrans, col=colSet02[binId2])
For this exercise only use the sp and raster packages (and indirectly, packages such as rgdal that raster may use in the background). Consider answering by writing sentences that use inline R commands like 5+5 = ‘r5+5’. a) Read the shapefile “galap.shp” (in galap.zip) into R to create a SpatialPolygonsDataFrame.
library(raster)
galap <- shapefile("galap.shp") #read galap data
class(galap)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
The class is a “SpatialPolygonsDataFrame”.
length(galap@polygons)
## [1] 30
There are 30 polygons in this object.
dim(galap)
## [1] 30 12
There are 12 variables.
plot(galap)
plot(galap$species~galap$area)
I would use species count divided by the log of area. Based on the scatterplot above, I wouldn’t want to simply divide number of species by units of area, because the data is not well distributed. Taking the log of the area normalizes the data, as shown in the scatterplot below.
plot(galap$species~log(galap$area))
galap$density = galap$species/log(galap$area)
spplot(galap, 'density')
which.max(galap$area) #query which polygon has the largest area
## [1] 2
galap$Island[2] #identify the name of this island
## [1] "Isabela"
Isabela=subset(galap, Island=='Isabela') #subset data to only include this island
plot(Isabela)
writeOGR(Isabela, datapath, 'Isabela_shp', driver='ESRI Shapefile') #write a new shapefile for this island
library(raster)
alt=getData('alt', country='Ecuador') #download altitude data for Ecuador
#change projection of galap to match that of alt
alt@crs #identify coordinate system of alt
## CRS arguments:
## +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
galap_longlat = spTransform(galap, crs(alt)) #transform galapagos shapefile
galap_crop = crop(alt, galap_longlat) #crop altitute to new galapagos shapefile
plot(galap_crop)
plot(galap_longlat, bg='transparent', add=T) #map elevation data with island outlines
contour(alt)